rm(list = ls(all = TRUE)) #clear workspace

setwd("~/Dropbox/REDISTRICTING/SIMULA_MVN")
source("02_get.parms.R")
source("03_simulate.R")
set.seed(12345) #Random seed for replication

# Computos 2012 con Redistritacion 2013
# Homologacion elaborada antes en la carpeta correspondiente
data <- read.table("DIPUTADOS_2012_2013.txt", header=TRUE, sep="\t")
data <- with(data, 
             data.frame(
               PAN=PAN,
               PRI=PRI+CPM,
               Mov.Prog=MP,
               OTH=PANAL + PVEM))

# Curules observadas
obs.seats <- apply(data, 1,function(x) max(x)==x)
obs.seats <- rowSums(obs.seats)

# Simula escenario hipotetico
#----------------------------
# Parametros
parms <- get.parms(data)

hip.mu <- c(.31,.31,.31,.07)
hip.mu <- log(hip.mu[-length(hip.mu)]/hip.mu[length(hip.mu)])
#hip.mu <- parms$mu
#hip.mu[1:3] <- mean(hip.mu[1:3])

#Variar mu, misma sigma ("swing ratio")
sim.hip <- replicate(1000, simulate(n=300, mu=hip.mu, Sigma=parms$sigma), simplify=FALSE)
sim.hip.seats <-do.call(rbind, sapply(sim.hip, "[", "seats"))

# Curules
apply(sim.hip.seats, 2, quantile, c(0.5,0.025,0.975))

